Cars arrive at the Costco gas station at rate λ = 1 car per minute.
- Cars arrive one-way to enter the gas station.
- There are 3 islands with 2 pumps on each side, for a total of 12 pumps. Label the pumps in the back from left to right as 1, 2, …, 6, and the pumps in front from left to right as 7, 8, …, 12.
- When cars arrive they join one of 6 queues, one queue for each side of the islands. There is one queue that waits for pumps 1 and 7, one that waits for pumps 2 and 8, and so on. The car first in line in the queue for pumps 1 and 7 will move to whichever one opens up first, similarly for 2 and 8, and so on.
- There are 3 different types of cars: 30% of cars will only join the queues on the left (to wait for one of the odd numbered pumps), 30% of cars will only join the queues on the right (to wait for an even numbered pump), and the remaining cars will join any queue. (Idea: the pumps are supposed to reach to either side of a car, but this is harder to do for some cars (e.g. minivans), so some cars will only join the queues corresponding to which side their gas tank is on.)
- When a car arrives, it will join whichever queue of its type is the shortest. If there is a tie for shortest queue of its type, the car will choose at random among the shortest queue options.
- Once a car chooses a queue, it does not switch to another or leave the gas station before completing service.
- Assume each pump serves customers at Exponential rate μ = 0.2 cars per minute, independent of the customer type. You can assume that once a car completes service the next car in the queue starts service immediately. (Essentially, you can assume that any time between when one car finishes pumping and the next one starts has been accounted for in the service rate.
Part 1
Creating Simulation
car <- setClass("car", slots = c(queue = "numeric",
start_time = "numeric", end_time = "numeric"))
setMethod("show", "car",
function(object) {
cat("Queue:", object@queue, "\n")
cat("Time in system:", object@start_time, "to", object@end_time)
})get_count <- function(lst, var, val) {
if(length(lst) == 0) {
return(0)
}
sum(sapply(lst, function(x) attr(x, var) == val))
}
filter_cars <- function(lst, var, val) {
lst[sapply(lst, function(x) attr(x, var) == val)]
}
sim <- function(t_end) {
W_n <- c(0)
X_t <- c(0)
lambda <- 1
mu <- 0.2
queues <- list()
num_queues <- rep(0, 6)
pumps <- list()
car_lst <- list()
while(W_n[length(W_n)] < t_end) {
# time next event happens
t <- W_n[length(W_n)] + rexp(1) / (lambda + mu * length(pumps))
# if next event happens after end time
if(t > t_end) {
W_n <- append(W_n, t_end)
X_t <- append(X_t, X_t[length(X_t)])
for(c in pumps) {
c@end_time <- -1
car_lst <- append(car_lst, c)
}
for(c in queues) {
c@end_time <- -1
car_lst <- append(car_lst, c)
}
return(list(W_n, X_t, car_lst))
}
W_n <- append(W_n, t)
# pick event
event <- sample(c(0, 1), 1, prob = c(mu * length(pumps), lambda))
# car arrives
if(event) {
# select queue
type <- sample(1:3, 1, prob = c(.3, .3, .4))
if(type == 1) {
min <- min(num_queues[1], num_queues[3], num_queues[5])
options <- sapply(c(num_queues[1], num_queues[3], num_queues[5]), function(x){x == min}) * c(1, 3, 5)
q <- sample(options[options != 0], 1)
} else if(type == 2) {
min <- min(num_queues[2], num_queues[4], num_queues[6])
options <- sapply(c(num_queues[2], num_queues[4], num_queues[6]), function(x){x == min}) * c(2, 4, 6)
q <- sample(options[options != 0], 1)
} else {
min <- min(num_queues)
options <- sapply(num_queues, function(x){x == min}) * 1:6
q <- sample(options[options != 0], 1)
}
# add to queue
queues <- append(queues, new("car", queue = q, start_time = t))
num_queues[q] <- num_queues[q] + 1
X_t <- append(X_t, X_t[length(X_t)] + 1)
# car leaves
} else {
# select pump
i <- sample(1:length(pumps), 1)
# remove from pump
pumps[[i]]@end_time <- t
q <- pumps[[i]]@queue
car_lst <- append(car_lst, pumps[[i]])
pumps <- pumps[-i]
X_t <- append(X_t, X_t[length(X_t)] - 1)
}
# add to pump if empty
if(num_queues[q] > 0 & get_count(pumps, "queue", q) < 2) {
pumps <- append(pumps, filter_cars(queues, "queue", q)[1])
i <- which(sapply(queues, function(x){x@start_time == pumps[length(pumps)][[1]]@start_time}))
queues <- queues[-i]
num_queues[q] <- num_queues[q] - 1
}
}
return(list(W_n, X_t, car_lst))
}results <- sim(10080)Results
Long run distribution of the number of cars in the system (waiting or in service)
df_results <- data.frame(time = results[[1]], count = results[[2]]) %>%
mutate(time_spent = lead(time) - time)
df_results[nrow(df_results), 3] = df_results[nrow(df_results), 1] - df_results[nrow(df_results) - 1, 1]
total_time_spent <- sum(df_results$time_spent)
df_results %>%
group_by(count) %>%
summarise(prop = sum(time_spent) / total_time_spent) %>%
ggplot() +
geom_col(aes(x = count, y = prop), color = "black", fill = "lightblue") +
theme_bw() +
labs(x = "Cars in System",
y = "Proportion of Time",
title = "Long Run Proportion of Time Spent with X Cars in System")
Long run fraction of time there are no cars in the system
kable(
df_results %>%
group_by(count) %>%
summarise(prop = sum(time_spent) / total_time_spent) %>%
filter(count == 0)
)| count | prop |
|---|---|
| 0 | 0.003361 |
Long run average number of customers in the system
mean(results[[2]])## [1] 6.029926Long run distribution of the amount of time a customer spends in the system (waiting time + service time)
df_cars <- data.frame(Queue = sapply(results[[3]], function(x){x@queue}),
Start = sapply(results[[3]], function(x){x@start_time}),
End = sapply(results[[3]], function(x){x@end_time})) %>%
filter(End != -1) %>%
mutate(Time = End - Start)
df_cars %>%
ggplot() +
geom_histogram(aes(x = Time, y = after_stat(density)), color = "black", fill = "lightblue") +
theme_bw() +
labs(x = "Time Spent in System",
y = "Proportion",
title = "Long Run Distribution of Time Spent in System")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Average time a customer spends in the system
avg_time <- mean(df_cars$Time)
avg_time## [1] 5.541578Proportion of cars in each queue
total_cars <- df_cars %>% count()
kable(
df_cars %>%
group_by(Queue) %>%
summarise(n() / total_cars)
)| Queue | n |
|---|---|
| 1 | 0.1679223 |
| 2 | 0.1663362 |
| 3 | 0.1666336 |
| 4 | 0.1670301 |
| 5 | 0.1628668 |
| 6 | 0.1692109 |
Bootstraped distribution of average number of cars in the system (just learned bootstrapping in 427!)
boot_avg <- numeric(10000)
for(i in 1:10000) {
boot_avg[i] <- mean(sample(results[[2]], length(results[[2]]), replace = TRUE))
}
hist(boot_avg)
Part 2
- Changes:
- If a car arrives and the shortest line is 10 cars or longer, it will leave.
- We have installed 3 super fast electric car charger (has its own queue). This charger serves customers at Exponential rate μ = 0.1 cars per minute. Electric cars will show up with probability 0.1 (any queue cars now 0.3). If the charger is occupied and an electric car shows up, it will leave.
sim2 <- function(t_end) {
W_n <- c(0)
X_t <- c(0)
lambda <- 1
mu <- 0.2
mu2 <- 0.1
queues <- list()
num_queues <- rep(0, 6)
pumps <- list()
charger <- list()
car_lst <- list()
while(W_n[length(W_n)] < t_end) {
# time next event happens
t <- W_n[length(W_n)] + rexp(1) / (lambda + mu * length(pumps) + mu2 * length(charger))
# if next event happens after end time
if(t > t_end) {
W_n <- append(W_n, t_end)
X_t <- append(X_t, X_t[length(X_t)])
for(c in pumps) {
c@end_time <- -1
car_lst <- append(car_lst, c)
}
for(c in charger) {
c@end_time <- -1
car_lst <- append(car_lst, c)
}
for(c in queues) {
c@end_time <- -1
car_lst <- append(car_lst, c)
}
return(list(W_n, X_t, car_lst))
}
# pick event
event <- sample(c(0, 1, 2), 1, prob = c(mu * length(pumps), lambda, mu2 * length(charger)))
# car arrives
if(event == 1) {
# select queue
type <- sample(1:4, 1, prob = c(.3, .3, .3, .1))
if(type == 1) {
min <- min(num_queues[1], num_queues[3], num_queues[5])
options <- sapply(c(num_queues[1], num_queues[3], num_queues[5]), function(x){x == min}) * c(1, 3, 5)
q <- sample(options[options != 0], 1)
} else if(type == 2) {
min <- min(num_queues[2], num_queues[4], num_queues[6])
options <- sapply(c(num_queues[2], num_queues[4], num_queues[6]), function(x){x == min}) * c(2, 4, 6)
q <- sample(options[options != 0], 1)
} else if(type == 3) {
min <- min(num_queues)
options <- sapply(num_queues, function(x){x == min}) * 1:6
q <- sample(options[options != 0], 1)
} else {
q <- 7
if(length(charger) < 3) {
charger <- append(charger, new("car", queue = 7, start_time = t))
W_n <- append(W_n, t)
X_t <- append(X_t, X_t[length(X_t)] + 1)
}
}
# add to queue
if(q != 7) {
if(num_queues[q] < 10) {
queues <- append(queues, new("car", queue = q, start_time = t))
num_queues[q] <- num_queues[q] + 1
W_n <- append(W_n, t)
X_t <- append(X_t, X_t[length(X_t)] + 1)
}
}
# gas car leaves
} else if (event == 0) {
W_n <- append(W_n, t)
# select pump
i <- sample(1:length(pumps), 1)
# remove from pump
pumps[[i]]@end_time <- t
q <- pumps[[i]]@queue
car_lst <- append(car_lst, pumps[[i]])
pumps <- pumps[-i]
X_t <- append(X_t, X_t[length(X_t)] - 1)
# electric car leaves
} else {
W_n <- append(W_n, t)
# select charger
i <- sample(1:length(charger), 1)
# remove from charger
charger[[i]]@end_time <- t
car_lst <- append(car_lst, charger[[i]])
charger <- charger[-i]
X_t <- append(X_t, X_t[length(X_t)] - 1)
}
# add to pump if empty
if(q != 7) {
if(num_queues[q] > 0 & get_count(pumps, "queue", q) < 2) {
pumps <- append(pumps, filter_cars(queues, "queue", q)[1])
i <- which(sapply(queues, function(x){x@start_time == pumps[length(pumps)][[1]]@start_time}))
queues <- queues[-i]
num_queues[q] <- num_queues[q] - 1
}
}
}
return(list(W_n, X_t, car_lst))
}results2 <- sim2(10080)Long run average number of customers in the system
mean(results2[[2]])## [1] 6.293181Average time a customer spends in the system
df_cars2 <- data.frame(Queue = sapply(results2[[3]], function(x){x@queue}),
Start = sapply(results2[[3]], function(x){x@start_time}),
End = sapply(results2[[3]], function(x){x@end_time})) %>%
filter(End != -1) %>%
mutate(Time = End - Start)
avg_time2 <- mean(df_cars2$Time)
avg_time2## [1] 5.812547Proportion of time each queue was in use by at least 1 car
df_results2 <- data.frame(time = results2[[1]], count = results2[[2]]) %>%
mutate(time_spent = lead(time) - time)
df_results2[nrow(df_results2), 3] = df_results2[nrow(df_results2), 1] - df_results2[nrow(df_results2) - 1, 1]
total_time_spent2 <- sum(df_results2$time_spent)
kable(df_cars2 %>%
group_by(Queue) %>%
summarise(sum(Time) / total_time_spent2)
)| Queue | sum(Time)/total_time_spent2 |
|---|---|
| 1 | 0.8134380 |
| 2 | 0.8419245 |
| 3 | 0.7919640 |
| 4 | 0.8053084 |
| 5 | 0.7860441 |
| 6 | 0.7814158 |
| 7 | 0.9356341 |
Long run fraction of time there are no cars in the system
kable(
df_results2 %>%
group_by(count) %>%
summarise(prop = sum(time_spent) / total_time_spent2) %>%
filter(count == 0)
)| count | prop |
|---|---|
| 0 | 0.0038369 |
- As expected, both averages increased due to the addition of 3 new “pumps.” It appears that cars leaving if the shortest queue is 10 cars or longer as well as not allowing queues for the chargers did not affect the averages as much as the addition of the new “pumps.” We can see that the chargers had a large impact on increasing the count of cars in the system since at least one charger was in use for 90% of the total time, thus raising the average amount of cars in the system.